//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Gauss_Legendre_1D.cc
 * \author Jeremy Roberts
 * \date   06/19/2011
 * \brief  Gauss_Legendre_1D member definitions.
 * \note   Copyright (C) 2011 Jeremy Roberts
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <cstdio>
#include <cmath>

#include "angle/Gauss_Legendre_1D.hh"

#include "utilities/DBC.hh"
#include "utilities/Soft_Equivalence.hh"
#include "utilities/Constants.hh"

namespace slabtran
{

//---------------------------------------------------------------------------//
// CONSTRUCTOR
//---------------------------------------------------------------------------//
/*!
 * \brief Constructor.
 */
Gauss_Legendre_1D::Gauss_Legendre_1D(int    sn_order,
		                             double norm)
: Base(sn_order, norm)
{
	using util::soft_equiv;

	Require (d_sn_order > 0);      // need a positive order
	Require (d_norm > 0.0);        // need a positive normalization
	Require (d_sn_order % 2 == 0); // need an even order

	// number of angles in the quadrature
	int numAngles = d_sn_order;

	// size angles (1-D)
	//d_angles.resize(1);

	// get references to mu
	Base::Vec_Dbl &mu  = d_angles;//[MU];

	// The number of quadrature levels is equal to the requested SN order.
	size_t levels       = d_sn_order;
	size_t octantAngles = numAngles / 2;  // 2 octants in 1D.

	// Force the direction vectors to be the correct length.
	mu.resize(numAngles);
	d_wt.resize(numAngles);

	// Set GL values for one half-space based on the SN order specified.
	// For 2:2:20, we use hard-coded, high precision values, though I honestly
	//  don't know the value in giving beyond 16 sig figs...anyone know?
	// For order > 20, we compute the parameters numerically
	switch (levels)
	{
		case 2:
		{
			// GL S-2 Quadrature Set.
			mu[0] = 0.5773502691896257645091488;
			d_wt[0] = 1.0000000000000000000000000;
			break;
		}

		case 4:
		{
			// GL S-4 Quadrature Set
			mu[0] = 0.8611363115940525752239465;
			mu[1] = 0.3399810435848562648026658;
			d_wt[0] = 0.3478548451374538573730639;
			d_wt[1] = 0.6521451548625461426269361;
			break;
		}

		case 6:
		{
			// GL S-6 Quadrature Set
			mu[0] = 0.9324695142031520278123016;
			mu[1] = 0.6612093864662645136613996;
			mu[2] = 0.2386191860831969086305017;
			d_wt[0] = 0.1713244923791703450402961;
			d_wt[1] = 0.3607615730481386075698335;
			d_wt[2] = 0.4679139345726910473898703;
			break;
		}

		case 8:
		{
			// GL S-8 Quadrature Set
			mu[0] = 0.9602898564975362316835609;
			mu[1] = 0.7966664774136267395915539;
			mu[2] = 0.5255324099163289858177390;
			mu[3] = 0.1834346424956498049394761;
			d_wt[0] = 0.1012285362903762591525314;
			d_wt[1] = 0.2223810344533744705443560;
			d_wt[2] = 0.3137066458778872873379622;
			d_wt[3] = 0.3626837833783619829651504;
			break;
		}

		case 10:
		{
			// GL S-10 Quadrature Set
			mu[0] = 0.9739065285171717200779640;
			mu[1] = 0.8650633666889845107320967;
			mu[2] = 0.6794095682990244062343274;
			mu[3] = 0.4333953941292471907992659;
			mu[4] = 0.1488743389816312108848260;
			d_wt[0] = 0.0666713443086881375935688;
			d_wt[1] = 0.1494513491505805931457763;
			d_wt[2] = 0.2190863625159820439955349;
			d_wt[3] = 0.2692667193099963550912269;
			d_wt[4] = 0.2955242247147528701738930;
			break;
		}

		case 12:
		{
			// GL S-12 Quadrature Set */
			mu[0] = 0.9815606342467192506905491;
			mu[1] = 0.9041172563704748566784659;
			mu[2] = 0.7699026741943046870368938;
			mu[3] = 0.5873179542866174472967024;
			mu[4] = 0.3678314989981801937526915;
			mu[5] = 0.1252334085114689154724414;
			d_wt[0] = 0.0471753363865118271946160;
			d_wt[1] = 0.1069393259953184309602547;
			d_wt[2] = 0.1600783285433462263346525;
			d_wt[3] = 0.2031674267230659217490645;
			d_wt[4] = 0.2334925365383548087608499;
			d_wt[5] = 0.2491470458134027850005624;
			break;
		}

		case 14:
		{
			// GL S-14 Quadrature Set
			mu[0] = 0.9862838086968123388415973;
			mu[1] = 0.9284348836635735173363911;
			mu[2] = 0.8272013150697649931897947;
			mu[3] = 0.6872929048116854701480198;
			mu[4] = 0.5152486363581540919652907;
			mu[5] = 0.3191123689278897604356718;
			mu[6] = 0.1080549487073436620662447;
			d_wt[0] = 0.0351194603317518630318329;
			d_wt[1] = 0.0801580871597602098056333;
			d_wt[2] = 0.1215185706879031846894148;
			d_wt[3] = 0.1572031671581935345696019;
			d_wt[4] = 0.1855383974779378137417166;
			d_wt[5] = 0.2051984637212956039659241;
			d_wt[6] = 0.2152638534631577901958764;
			break;
		}

		case 16:
		{
			// GL S-16 Quadrature Set
			mu[0] = 0.9894009349916499325961542;
			mu[1] = 0.9445750230732325760779884;
			mu[2] = 0.8656312023878317438804679;
			mu[3] = 0.7554044083550030338951012;
			mu[4] = 0.6178762444026437484466718;
			mu[5] = 0.4580167776572273863424194;
			mu[6] = 0.2816035507792589132304605;
			mu[7] = 0.0950125098376374401853193;
			d_wt[0] = 0.0271524594117540948517806;
			d_wt[1] = 0.0622535239386478928628438;
			d_wt[2] = 0.0951585116824927848099251;
			d_wt[3] = 0.1246289712555338720524763;
			d_wt[4] = 0.1495959888165767320815017;
			d_wt[5] = 0.1691565193950025381893121;
			d_wt[6] = 0.1826034150449235888667637;
			d_wt[7] = 0.1894506104550684962853967;
			break;
		}

		case 18:
		{
			// GL S-18 Quadrature Set
			mu[0] = 0.9915651684209309467300160;
			mu[1] = 0.9558239495713977551811959;
			mu[2] = 0.8926024664975557392060606;
			mu[3] = 0.8037049589725231156824175;
			mu[4] = 0.6916870430603532078748911;
			mu[5] = 0.5597708310739475346078715;
			mu[6] = 0.4117511614628426460359318;
			mu[7] = 0.2518862256915055095889729;
			mu[8] = 0.0847750130417353012422619;
			d_wt[0] = 0.0216160135264833103133427;
			d_wt[1] = 0.0497145488949697964533349;
			d_wt[2] = 0.0764257302548890565291297;
			d_wt[3] = 0.1009420441062871655628140;
			d_wt[4] = 0.1225552067114784601845191;
			d_wt[5] = 0.1406429146706506512047313;
			d_wt[6] = 0.1546846751262652449254180;
			d_wt[7] = 0.1642764837458327229860538;
			d_wt[8] = 0.1691423829631435918406565;
			break;
		}

		case 20:
		{
			// GL S-20 Quadrature Set
			mu[0] = 0.9931285991850949247861224;
			mu[1] = 0.9639719272779137912676661;
			mu[2] = 0.9122344282513259058677524;
			mu[3] = 0.8391169718222188233945291;
			mu[4] = 0.7463319064601507926143051;
			mu[5] = 0.6360536807265150254528367;
			mu[6] = 0.5108670019508270980043641;
			mu[7] = 0.3737060887154195606725482;
			mu[8] = 0.2277858511416450780804962;
			mu[9] = 0.0765265211334973337546404;
			d_wt[0] = 0.0176140071391521183118620;
			d_wt[1] = 0.0406014298003869413310400;
			d_wt[2] = 0.0626720483341090635695065;
			d_wt[3] = 0.0832767415767047487247581;
			d_wt[4] = 0.1019301198172404350367501;
			d_wt[5] = 0.1181945319615184173123774;
			d_wt[6] = 0.1316886384491766268984945;
			d_wt[7] = 0.1420961093183820513292983;
			d_wt[8] = 0.1491729864726037467878287;
			d_wt[9] = 0.1527533871307258506980843;
			break;
		}

		default:
		{
			// generate the parameters on-the-fly
			double tmpmu[numAngles];
			double tmpwt[numAngles];
			Generate_Parameters(numAngles,tmpmu,tmpwt);
			for ( int i = 0; i < numAngles/2; ++i )
			{
				mu[i]   = tmpmu[i];
				d_wt[i] = tmpwt[i];
			}
			break;
		}
	}

	// Evaluate mu and wt for second half-space
	for ( int i = 0; i < numAngles/2; ++i )
	{
		mu[numAngles/2+i]   = -mu[numAngles/2-i-1];
		d_wt[numAngles/2+i] = d_wt[numAngles/2-i-1];
	}

	// Normalize the quadrature set
	double wsum = 0.0;
	for (size_t n = 0; n <= numAngles - 1; ++n)
		wsum = wsum + d_wt[n];

	for (size_t n = 0; n <= numAngles - 1; ++n)
		d_wt[n] = d_wt[n] * (norm / wsum);


#ifdef DEBUG
	display();
#endif

} // end constructor

/*!
 * \brief Generate Gauss-Legendre parameters.
 *
 * Modified version of the function gauleg from <B> Numerical Recipes in C</B>
 * (Cambridge Univ. Press) by W.H. Press, S.A. Teukolsky, W.T. Vetterling, and
 * & B.P. Flannery
 *
 * \param n  number of points (i.e. the quadrature order)
 * \param x  temporary array for Gauss points
 * \param w  temporary array for weights
 *
 */
void Gauss_Legendre_1D::Generate_Parameters( int n, double *x, double *w )
{
	int m,j,i;
	double z1,z,xm,xl,pp,p3,p2,p1;
	m=(n+1)/2; /* The roots are symmetric, so we only find half of them. */
	for (i=1;i<=m;i++) { /* Loop over the desired roots. */
		z=cos(constants::pi*(i-0.25)/(n+0.5));
		/* Starting with the above approximation to the ith root, we enter */
		/* the main loop of refinement by Newton's method.                 */
		int count = 0;
		do {
			p1=1.0;
			p2=0.0;
			for (j=1;j<=n;j++) { /* Recurrence to get Legendre polynomial. */
				p3=p2;
				p2=p1;
				p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
			}
			/* p1 is now the desired Legendre polynomial. We next compute */
			/* pp, its derivative, by a standard relation involving also  */
			/* p2, the polynomial of one lower order.                     */
			pp=n*(z*p1-p2)/(z*z-1.0);
			z1=z;
			z=z1-p1/pp; /* Newton's method. */
		} while ( std::abs(z-z1) > 1e-15 && ++count < 200 );
		x[i-1]=z;                     /* Store root */
		w[i-1]=2.0/((1.0-z*z)*pp*pp); /* Compute the weight */
	}
}


//---------------------------------------------------------------------------//
// PUBLIC FUNCTIONS
//---------------------------------------------------------------------------//
/*!
 * \brief Display the quadrature.
 */
void Gauss_Legendre_1D::display() const
{

	using std::cout;
	using std::endl;
	using std::printf;

	cout << endl << "The Quadrature directions and weights are:"
			<< endl << endl;
	cout << "   m         mu              wt       " << endl;
	cout << "  ---   --------------  --------------" << endl;
	double sum_wt = 0.0;
	for ( int i = 0; i < d_angles.size(); ++i )
	{
		printf ("%4i    %13.10f   %13.10f \n", i, d_angles[i], d_wt[i]);
		sum_wt += d_wt[i];
	}
	cout << endl << "  The sum of the weights is " << sum_wt << endl;
	cout << endl;
}

} // end namespace slabtran

//---------------------------------------------------------------------------//
//                 end of Gauss_Legendre_1D.cc
//---------------------------------------------------------------------------//

